function outdat = fsb_detrend_nan(indat,poly_order,gopt)

% FSB: remove linear trends from data that contains NaN values
%
% EXAMPLE:
% idat = fsb_detrend_nan(idat,5,0,0)
%
% INPUT:
% indat: data to be processed
% poly_order: polygon order
% gopt: display graphical output
%
% OUTPUT:
% idat: detrended 4-D image data
%
% CALLED BY:
% fsb_detrend_scan.m
%
% NOTES:
% Accepts vectors and 2-D arrays
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
%
%$Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<3;
    gopt =1;
end

dim2 = min(size(indat));
dim = size(indat);

switch dim2

    case 1 % Data is provided as a vector

        indat = double(indat(:));
        y = [1: size(indat,1)];

        %~~~~~~~~~~~~~~~~~~~
        % Find NaN values
        %~~~~~~~~~~~~~~~~~~~
        nanny = find(isnan(indat));
        indat(nanny) = [];
        y(nanny) = [];

        %~~~~~~~~~~~~~~~~~~~
        % Calculate fit parameters
        %~~~~~~~~~~~~~~~~~~~
        [p,ErrorEst] = polyfit(y,indat',poly_order);

        %~~~~~~~~~~~~~~~~~~~
        % Evaluate the fit
        %~~~~~~~~~~~~~~~~~~~
        corrvec1 = polyval(p,y,ErrorEst);
        corrvec2 = mean(indat(:))-corrvec1;
        outdat1 = indat'+corrvec2;
        [p,ErrorEst] = polyfit(y,outdat1,poly_order);
        endcorr = polyval(p,y,ErrorEst);

        %~~~~~~~~~~~~~~~~~~~
        % Put it back into right place
        %~~~~~~~~~~~~~~~~~~~
        outdat = nan(dim);
        outdat(y) = outdat1;
        indat2 = nan(dim);
        indat2(y) = indat;

        if gopt ==1;

            %~~~~~~~~~~~~~~~~~~~
            % Plot the data and the fit
            %~~~~~~~~~~~~~~~~~~~
            figure(1); cla;
            subplot(2,1,1)
            plot(y,corrvec1,'-',y,indat','+',y,outdat1','r.',y,endcorr,'m');
            legend('Polynomial Model','Original Data','Corrected data','Fit after correction');
            xlabel('Volume');
            ylabel('Signal intensity');

            %~~~~~~~~~~~~~~~~~~~
            % Plot the original and detrended data
            %~~~~~~~~~~~~~~~~~~~
            subplot(2,1,2);
            y = [1: size(indat2,1)];
            plot(y,outdat,'r',y,indat2,'g');
            legend('Corrected data','Original Data');
            xlabel('Volume');
            ylabel('Signal intensity');
        end

    otherwise %Data is provided as an array

        tic;
        indat2d = single(indat);
        outdat = nan(dim);
        Y = 1: size(indat,1);
        h = waitbar(0,'Detrending data');

        for x = 1: size(indat2d,2)
            waitstring = ['Detrending data : ' num2str(int16(toc)) ' seconds elapsed'];
            waitbar(x/size(indat2d,2),h,waitstring);

            %~~~~~~~~~~~~~~~~~~~
            % Take indat vector
            %~~~~~~~~~~~~~~~~~~~
            indat = indat2d(:,x);
            y = 1: size(indat,1);

            %~~~~~~~~~~~~~~~~~~~
            % Find NaN values
            %~~~~~~~~~~~~~~~~~~~
            nanny = find(isnan(indat));
            indat(nanny) = [];
            y(nanny) = [];

            %~~~~~~~~~~~~~~~~~~~
            % Calculate fit parameters
            %~~~~~~~~~~~~~~~~~~~
            [p,ErrorEst] = polyfit(y,indat',poly_order);
            % Evaluate the fit
            corrvec1 = polyval(p,y,ErrorEst);
            corrvec2 = mean(indat(:))-corrvec1;
            outdat1 = indat'+corrvec2;

            %~~~~~~~~~~~~~~~~~~~
            % Put it back into right place
            %~~~~~~~~~~~~~~~~~~~
            outdat(y,x) = outdat1';

        end
        close(h);
        outmean = mean(outdat,2);
        inmean = mean(indat2d,2);
        t=toc; disp(['Time elapsed ' num2str(t) ' s']);

        if gopt ==1
            %~~~~~~~~~~~~~~~~~~~
            % Plot the data and the corrected data
            %~~~~~~~~~~~~~~~~~~~
            figure(1); clf;
            plot(Y,inmean,'.r',Y,outmean','+');
            legend('Original Data mean','Corrected data mean');
            xlabel('Volume');
            ylabel('Signal intensity');
        end

end

end
